A model for the atomic-scale structure of a dense, nonequilibrium fluid: the 
homogeneous cooling state of granular fluids 
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It is shown that the equilibrium Generalized Mean Spherical Model of fluid structure may be 
(3JT)' extended to nonequilibrium states with equation of state information used in equilibrium replaced by 

an exact condition on the two-body distribution function. The model is applied to the homogeneous 
cooling state of granular fluids and upon comparison to molecular dynamics simulations is found to 
provide an accurate picture of the pair distribution function. 
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^ ' I. INTRODUCTION 

The atomic scale structure, e.g. the pair distribution function or, equivalently, the density-density equal time 
correlation function, of both equilibrium and nonequilibrium fluids is directly accessible to experiment by means of 
' light scattering and has been used to study the behavior of complex systems such as sheared colloidal suspensions 
jij. However, the subject of nonequilibrium fluid structure remains obscure, particularly when compared to the 
structure of equilibrium fluids which is one of the most advanced areas of equilibrium statistical mechanics [|| . While 
kinetic theory can provide some asymptotic results, see e.g. ||, 0], models of the small-scale structure remain 
phenomenological. This is not surprising: in equilibrium, one starts with the exact N-body distribution and the 
| problem is to integrate out uninteresting degrees of freedom so as to get the pair distribution function (pdf). Away 
i from equilibrium, even in steady states, the N-body distribution function is not known and the only starting point is 
either a complex dynamical equation, the Liouville equation, or equivalently a coupled set of equations, the BBGKY 
hierarchy relating the n-body to the (n+l)-body distribution functions. Another method used to derive equilibrium 
integral equations is based on the fact that from the BBGKY hierarchy it is easy, essentially by integrating over 
the velocities, to derive the YBG hierarchy which relates the pair distribution function to the triplet distribution 
function, the triplet to the quartic distribution function, etc. Some sort of closure hypothesis, like the neglect of three 
body correlations, then yields a closed equation for the pdf, see e.g. [||, ||. However, away from equilibrium, the 
"p*; ■ distribution of velocities is not known and this procedure cannot be followed in detail. One approach to nonequilibrium 
fluid structure is that of Hess and coworkers jj, ||. The object there was to describe sheared colloidal suspensions 
as well as computer simulations of sheared simple fluids. The theory, based on a kind of relaxation approximation, is 
phenomenological and will be considered in more detail below. Fluctuating hydrodynamics offers another possibility 
||, fLOfl , but is restricted to length-scales much larger than atomic length scales. Another approach, developed 
about a decade later by Eu et al |l2| , [|l3| , was to try to generalize the equilibrium integral equations, like the Kirkwood 
and the Percus-Yevick equations , to nonequilibrium fluids. For the reasons just mentioned, the only way to carry out 
this program was to start with an ansatz for the n-body distribution which allows the development to parallel that 
of the equilibrium equations. Whether or not this is a good approach is difficult to assess since the approximations 
made are uncontrolled, but it has been previously shown fl4|| , and will be argued in more detail below, that velocity 
correlations are the crucial determinant of atomic-scale deviations from equilibrium structure and, it would follow, 
that approaches which do consider them may be missing an essential contribution to the structure. Indeed, the role 
played by velocity correlations in determining one aspect of the structure of sheared simple fluids, namely the value 



and angular dependence of the pair distribution function at contact has already been demonstrated |15| and the 
purpose of this paper is to show that the knowledge of such velocity correlations can be used to extend equilibrium 
structural models into the nonequilibrium domain. 

In order to avoid a number of complexities associated with the spatial anisotropy of sheared fluids, the system used 
here to illustrate this model is the dissipative hard sphere model used in the study of granular fluids. This fundamen- 
tally nonequilibrium system presents a unique opportunity for the application and development of nonequilibrium 
statistical mechanics in a system which is both of practical relevance and yet sufficiently simple to be amenable to 
theoretical analysis. The basis for such analysis has recently been formulated by Brey, Dufty and Santos jl6| who have 
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constructed a Liouville description of the dynamics of the dissipative hard-sphere model. From this, Boltzmann- and 
Enskog-level kinetic equations for the one-body distribution function (the dense-fluid generalization of the Boltzmann 
equation) follow and their properties have been studied theoretically, in particular, the Chapman-Enskog solution to 
it has been developed, giving explicit expressions for the pressure and transport properties of the system [pL 7|| , and it is 
possible to solve it numerically p8[ . Despite the simplicity of the dynamics, simulations have shown that this system 
exhibits a rich phenomenology which is only partially understood at present (see, e.g., |L9|| ). The starting point for 
describing the phenomenology is the so-called homogeneous cooling state (HCS): since the collisions dissipate energy, 
the analog of the equilibrium state is one which is spatially homogeneous but in which the temperature decreases 
algebraically with time. However, this state is relatively unstable as it is subject to a number of forms of spontaneous 
symmetry breaking. The simplest example of this is the so-called clustering instability in which sufficiently large 
systems, subject to no external forces, exhibit a spatial clustering. This has been shown to be due to a small wave- 
vector hydrodynamic instability ]20| . Small systems for which the maximum wavevector (as determined by the size 
of the simulation cell) is too large, do not exhibit the transition. Other phenomena include the "shearing state", in 
which the system remains structurally homogeneous but long-ranged momentum correlations develop, and "inelastic 
collapse" in which a few particles collide with one-another many times until all relative momentum is dissipated and 
the particles come to rest in contact with one another. One motivation for trying to understand the structure of the 
HCS is that it is needed as input to various kinetic theory calculations which could be useful in understanding some 
of these features. 

The remainder of this paper is organized as follows. In Section II, the relevant elements of kinetic theory are reviewed 
and it is shown that the pdf at contact can in general be explicitly calculated with the same level of approximation as 
is used to get the Enskog equation. It is noted that both the pdf at contact, as well as the pressure, can be calculated 
in the HCS with no further approximations even though the exact solution for the one-body distribution is not known. 
In Section III, the dynamics of the dissipative hard sphere model is reformulated, through a change of variables, so 
that the HCS is mapped onto a steady state. Although mathematically, equivalent, the steady state formulation is 
convenient since it avoids numerical problems associated with the rapid decay of kinetic energy in the HCS and also 
allows statistical properties, like the pressure, to be computed by replacing ensemble averages by time averages (the 
assumption of ergodicity). Section IV consists of a presentation of a model for the structure of the HCS and Section 
V presents a comparison of this model to the results of simulations. The paper concludes with a discussion of the 
results and their bearing on the questions raised above. 



II. THEORY 



The model granular fluid consists of a collection of TV classical hard-spheres of diameter a having positions q and 
momenta p which experience inelastic binary collisions. Between collisions, the atoms freely stream so the state 
evolves according to: 



mi—qi = pi 
at 

d 

dt Pi = Fi 



(1) 
(2) 



where we allow for the possibility of an external force acting on the particles. Henceforth, we will consider only 
identical atoms and will use units in which the mass is equal to one. When two particles, say i and j, collide, their 
total momentum is unaffected but their relative momentum is altered according to 



Pij > hjPij = Pij - (1 + a) pij ■ % 



(3) 



where, in general, relative quantities are denoted as py = p.; — Pj , unit vectors as q = q/ |q| and where this equation 
serves to define the momentum transfer operator fey . The parameter a is called the coefficient of restitution and 
takes on values between one (elastic hard spheres) and zero (completely dissipative inelastic collisions). Under this 
dynamics, the one-body distribution function, fi(x, t) = /i (q, p, t) and the two-body distribution function ji(x\,xi,t) 
are related by the first equation of the BBGKY hierarchy Jlq] 



d . , . d „ . . 

g^fi[Xi,t) +pi • — h[xi,t) 



dxi I da S (qi2 - aa) (a ■ p i2 ) 



9 (-qi2 • P12) h{xi,x 2 ,t) 



(4) 
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where the notation indicates that integration of a is over the unit sphere, Q(x) — 1 for x > and is otherwise zero 
and where is the inverse of bij and may be easily seen to be 

Kj~Pi 3 = Pij - (1 + aT 1 ) Pij • (5) 



Finally, the two-body distribution must satisfy the identity, a kind of boundary condition 15 1, 

5 (qi2 - ad) 9 (qi2 • Pi 2 ) h{x\,xi,t) = 8 (q i2 - ad) 6 (q i2 • P12) -sb^/afci, x 2 , t) 

a 

= <J(qi2 - c?) -jb^O (-qi2 -Pi 2 )/ 2 (2;i,a; 2 ,t) (6) 
a 

the origin of which - basically, the conservation of probability during a collision - is discussed in the appendix. Here, 
I only note that this identity is independent of, and provides information additional to, the first BBGKY hierarchy 
equation given above. In fact, as discussed in the appendix, it can actually be derived from the second equation of 
the BBGKY hierarchy. Using an abbreviated notation, W\ 2 = 5 (q 12 — era), it is easily seen that the full distribution 
at contact can be written as 

Wuf2(xi,X2,t) = M / i 2 9(-qi 2 • p 12 ) f 2 (xi,x 2 ,t) +Wi 2 6(qi2 • P12) h{xi,x 2 ,t) 

= W 12 <d (-qi2 • P12) h{xi 1 x 2 ,t) + W 12 \b^ 2 Q (-q X2 • P12) h{xi,x 2 ,t) (7) 

or 

The combination 5 (qi 2 — 00) (— qi 2 • P12) f 2(^-1, ^2, t) appearing here, as well as in the collisional term in eq.(^) 
refers to the probability for two atoms just prior to a collision and we refer to it as the pre-collisional part of the 
distribution. Equation (|J) thus expresses the two-body distribution at contact solely in terms of the pre-collisional 
distribution. As discussed in ref. [fl5j| , the assumption of "molecular chaos" used to obtain the Boltzmann equation 
and its generalization to dense hard-sphere fluids, the Enskog equation, is that this pre-collisional distribution can be 
approximated by neglecting velocity correlations between the particles so that one writes 

W 12 Q (-qi2 • P12) f2(xi,x 2 ,t) ~ W 12 Q (-qia • P12) flo(qi, q2)/i(aci, t)fi(x 2 , t) (8) 

where go(qi 5 q 2 i^) is the pre-collisional pair distribution function (pdf) at time t which is normally taken to be the 
local-equilibrium form [pl| , E^j. Using this in eq. (Q), the generalized Enskog equation immediately results 

d . . , d . . . 

Q^h{xi,t) +pi • — }i{xi,t) 

= a 2 J dx 2 J da S (q i2 - ad) (a ■ p 12 ) a" 2 ^ 1 + 1 9 (-qi 2 ■ pi 2 ) 5o(qi, q 2 )/i(xi, t)f 1 (x 2 , t) (9) 

while from the boundary condition, we get 

Wi 2 f 2 {x 1 ,x 2 ,t) = W 12 Q (-q 12 • p 12 ) go(qi,<l2)fi(xi,t)fi{x2,t) 

+W 12 —b^ 2 1 Q (-qi 2 • pi 2 ) 5o (qi,q 2 )A(a;i,t)/i(a; 2 ,t) 
= M /r i 2 .9o(qi,q 2 )/i(^i^)/i(^ 2 ,i) 

+W 12 (q 12 • p 12 ) (Jpb£ - lj <7o(qi, q*)fi(*i> *)/i(^2, t) (10) 

thus expressing the distribution at contact as the sum of an uncorrelated term, the first on the right, and a term 
expressing the corrections due to momentum correlations. In equilibrium, the second term vanishes. Naturally, the 
momenta of the atoms are correlated after a collision and, in fact, ( ^tj| ) allows us to determine these correlations in a 
manner consistent with the degree of approximation of the generalized Enskog equation. Although such an evaluation 
might be used in a number of different ways, I will here focus one particular application of it to the problem of 
understanding the structure of the non-equilibrium state. Specifically, integrating over momenta gives 

Wi 2 n(qi; i)ra(q 2 ; t)g(qx, q 2 ; t) 
= IUi 2 n(qi)n(q 2 )5o(qi,q 2 ;£) 

+tUi 2 5o(qi,q 2 ) Jdp 1 dp2&{qx2 -Pis) (J^^ - A fi(xi,t)fi(x 2 ,t) (11) 
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where the nonequilibrium density is n(qi;i) = J dpifi(xi,t). Eq. ( |Tl| ) thus gives us an approximate evaluation of 
the contribution of momentum correlations to the structure of the fluid as expressed through the pdf. This relation 
was first used in ref. |ll| to characterize velocity correlations in a sheared fluid and has recently been used in a study 
of velocity correlations in granular fluids near equilibrium [ ^3| . 

The macroscopic balance equations for the density, momentum field mU, and energy density ^k B T follow from (^) 
by multiplying by 1, pi, and -^p\ respectively and integrating over p! to get (L7J 



d 

—n + V • nU = 
at 

£u t + U • VUi + (mn)" 1 ^ = 



dt 



^-T + U ■ VT + -L-pydjUi + V • q) = -T 3 / 2 C (12) 
at inks 

and explicit forms of the pressure tensor and heat flux vector q are given in the literature . The source term 
in the temperature equation is due to the cooling caused by the inelastic collisions. It is easy to see that a spatially 
homogeneous solution to these equations is possible with 

n(r, t) = n 
U = 

d_ T = _ T 3/2 C (13) 



where no is a constant. Because there is no potential energy in hard-sphere systems, the time and energy scales are set 
by the temperature and so, on dimensional grounds, it is clear that ( is independent of temperature: the temperature 

1 1/2 2 

is therefore given by T(t) = T 1 + iC^o * thus explicitly demonstrating the cooling. A linear stability analysis 

of the equations ( |l2|) expanded about the homogeneous cooling state shows that the state is unstable against small 
wavevector fluctuations p(| . 

The pdf at contact for HCS can also be evaluated without explicit knowledge of the one body distribution (see the 
Appendix) and is 



X(^)=yj dq_ 1 dq 2 S(q 1 2-aa)g(q 1 ,q2;t) 



= (M) 

where xo is defined by an analogous expression to this with ,g(qi, q2! t) replaced by go(qij c l2]t). Finally, since it will 
be of use below, we quote the expression for the collisional contribution to the pressure which is 



■ B T F 3 ^ nk B T 13 

i — 1 

= — j— 0-3 J da J dxidx 2 S(q 12 - a) 6 (q 12 • p i2 ) (qi2 • P12) 2 fi(xi, t)fi(x 2 , t)g (qx, c^) 

1 + a . . 

= ^^XO (15) 



where the last line follows from the spatial isotropy and time invariance of the HCS and also does not require explicit 
knowledge of the one-body distribution. This is useful since an exact solution to the Enskog equation for the one-body 
distribution for HCS is not known. 



III. MAPPING HCS TO A STEADY STATE 



Since the time scale is determined by the temperature, we can simplify the description of this state by making 
a change of variables. Specifically, define a scaled time coordinate via ws = ln(t/to) where w and to are arbitrary 
constants, and the corresponding equations of motion are 
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d_ 

ds 
d_ 

ds 



(16) 



where Ci — wtoe ws pi is the momentum (velocity) in the new coordinates. To carry through the statistical description, 
we must define a new set of probability densities. If we denote the set of positions and velocities of atoms x\...x m 
in the original system by T m and of those in the scaled system by T' m then the m-body distribution will transform 
according to 



(17) 



where J = 



9T' 



(wtoe ws ) mD is the Jacobian of the transformation for a D-dimensional system. We then find, 



e.g., that eq.(Q) becomes 
d 



~ d ~ d ~ 

g s H ( x l , S) + Cl • -Q^-Jl ( x l , s ) + ■ WCif! (x[ , s) 



a 2 / dx\dx' 2 / da S (qi2 — aa) (a ■ C12) 



1 



6(-qi2 ■ c 12 ) f 2 {x' l x' 2 ,t) 



(18) 



so that the only change is that a new term appears in the streaming operator. This term has the same form as the 
artificial thermostats used in the study of shear flow [jioj: however, here it arises solely from a change of variables. 
The balance equations become 

d „ 

— n + V • nV = 

ot 

^ Vi + V • W, + [mn)- x djPn - wVi = 

^-T + V • VT + ^-(PijdjVi + V • q) - 2wT = -T 3 / 2 ( (19) 
ot inks 

where V(q, s) = J dc c/ 1 (x' 1 ,s), |/csT(q, s) = J dc ^mc 2 fi(x[, s), etc. Now, the temperature of the homogeneous 
state is given by 




-1 e" 



(20) 



so that any initial temperature will equilibrate to a final stable value of {^^j ■ Thus, the original homogeneous 

cooling state is mapped by this change of variables onto a superficially steady state. It is worth noting that if one 
were to return to the original form of the dynamics and were to periodically rescale the momenta of the atoms so as 
to restore the temperature to its original value, then it is easy to see that as the time interval between rescalings goes 
to zero, the resulting equations of motion can be written in the form [l6| with the constant replaced by a complicated 
function of the momenta. Periodic rescaling is commonly used in simulations of sheared fluids and has recently been 
employed in the simulation of HCS p^l , 23 1. 

This mapping also makes the hydrodynamic instability apparent. If we expand (|l9| ) about the steady state and 
retain terms to second order in the gradients (i.e., the Navier-Stokes equations) and to first order in the densities and 
transform to Fourier space (with wavevector k), it is easy to see that the vorticity, w = kxV, satisfies 



d t ui + (Vofc 2 — w) ui = 



(21) 



where vq is the shear viscosity evaluated at density no an d temperature To = ( ^pj • It is obviously unstable for 
sufficiently small wavevectors and for sufficiently large systems, we therefore expect a spontaneous shear to develop 

(note that vq ~ \J~To so that the arbitrary constant, u>, plays no role in the stability criterion). It should also be 
noted that there is a closely related instability in the total (k = 0) velocity in the mapped system which, in the linear 
stability analysis, obeys 
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^Vi(p)-wVi(O) = (22) 

and so is clearly unstable for all system sizes. Since we are at liberty to choose the initial conditions to be those 
for which Vi(0) = 0, we will find that this is represents a minor problem in the simulations and is of no theoretical 
significance. 

IV. MODELLING THE STRUCTURE OF THE HCS 

In this section, we will assume spatial homogeneity so that the pdf depends only on the scalar separation between 
atoms. The simplest realistic model for the structure of the equilibrium hard sphere fluid is the Percus-Yevick 
approximation(see e.g., refs. p5|). This consists of the Ornstein-Zernike equation 

h(r) = c(r) + p { dr' c(|r - r'\)h(r') (23) 



where h(r) = g(r) — 1 and c(r) is the direct correlation function, together with the boundary conditions 

0(<7 - r)h{r) = -1 (24) 
Q(r-a)c(r) = 

where the first condition is exact while the second defines the approximation. Comparison with computer simulation 
shows that the PY approximation for the pair distribution function is quite accurate for separations greater that 
about two hard-sphere diameters but less accurate near contact. The description of the small-separation structure 
can be significantly improved by considering the Yukawa closure for the Orstein-Zernicke equation which replaces the 
boundary condition on the direct correlation function by 

i=l 

and choosing the constants K t and i?j to reproduce known properties: for example, taking m = 1 and fitting the 
Carnahan-Starling equation of state as calculated by both the pressure equation and the compressibility equation. 
The original Mean Spherical Approximation, for arbitrary pair potentials <&(r), consists of requiring that (d(r-cr)c(r) = 
<£>(r) where the effective hard sphere diameter is fit according to some criterion. The PY approximation is then seen 
to be the MSA for hard spheres. Equation (^||) may therefore be viewed as either the MSA for a potential which is 
the sum of Yukawas or as a general expansion with coefficients to be fitted in which case it is termed the Generalized 
MSA or GMSA and can be viewed as being systematic since any function could be fitted as a sum of Yukawas. 

Note that eq. ( p3| ) defines the direct correlation function and that the first of the boundary conditions in eq. (|24|) is 
an exact requirement. The only way in which this model uses the assumption that the system being modeled is in 
equilibrium is through the arguments that lead to the conclusion that the direct correlation function is short-ranged 
and hence the justification for the boundary conditions in eq.(p4]) and eq.(p5|). This connection to the equilibrium 
state is made even weaker in a reformulation of the model due to Yuste and Santos pq ], p7[ ], |2S[ ]. They begin by 
noting that the Laplace transform of the quantity rg(r), in the PY approximation, is naturally written as 



(26) 



G{t) = / dr e- sr rg{r) 



tF(t) 



-i 



1 - 127?F(£)e-* 
with 

m = ^Ts^Ts^Ts^ (27) 

They go on to point out that given the second equality of eq. (p6|) and making a Pade' approximation for the function 
F(t) one can deduce the correct order of the numerator and denominator of F(t) as well as the PY expressions for the 
coefficients based solely on the asymptotic properties of the pdf. Specifically, they note that a)g(r) at contact is given 
by 5( cr ) — linit_>oo t 2 F(t) thus fixing the relative number of terms in the numerator and denominator; b) linv^oo g(r) = 



G 



1 implies that G(t) — » t 2 and c) the fact that the static structure factor, given by S(q) = \im t ^i q Re (tG(t)), is finite 



t-> o 

at q = implies, given the previous condition, that for small t, G(t) — t~ 2 + o(l). (The last condition is equivalent 
to assuming that long-range correlations do not exist.) The minimal approximant satisfying these conditions is that 
given in eq.(|27|) with the PY value for the coefficients. They also note that the extension of the Pade' approximant to 
include one more term in both the numerator and denominator is exactly equivalent to the one- Yukawa closure while 
the further extension of the approximation corresponds to a closure consisting of a sum of Yukawas. This method is 
also shown to give the PY solution for sticky hard spheres as well as the exact structure for both ordinary and sticky 
hard spheres in one-dimension. A straight forward extension of these ideas has also been used to model the square-well 
fluid. Thus, in this formulation, the PY form of G(t) is taken as an ansatz characteristic of hard-core systems and 
the function F{t) modeled as a Pade' approximant subject to whatever knowledge exists about the structure. 

With this justification in mind, we consider the application of this approach to the nonequilibrium HCS. In equi- 
librium, the next inclusion of an additional term in the numerator and denominator of F(z) introduces two new 
parameters which are used to fit a known equation of state (normally the Carnahan-Starling equation of state) 
through both the pressure equation and the compressibility equation. The pressure equation for hard spheres in three 
dimensions reads 

P = 1 + ^VXeq (28) 



and so allows to calculate the pdf at contact, recall Xeq = 3(f), from the equation of state. In a nonequilibrium 
state, the collisional boundary condition can be used, together with assumption of molecular chaos, to give the same 
information. The equilibrium compressibility equation is 

|V) =l + pjdr (g(r) - 1) 

= l + km(G(t)-i) (29) 

and for this, there is no obvious substitute for the nonequilibrium state. With nothing to use in its place, I will 
continue to apply this even in the nonequilibrium state, calculating the pressure from eq. (|l5|) , which might be viewed 
as a local-equilibrium approximation. In fact, the resulting model is relatively insensitive to the value used for the 
pressure since this only fixes the area of the structure function whereas the results are quite sensitive to the value of 
the pdf at contact. Using go(qi,q2;i) ~ 9eq{q_i2) as is normally done in Enskog theory, the model is given by 

Fm _ 1 i + A 1 t + A 2 e 

W 127? So + Sit + S 2 t 2 + S 3 P + 1 ' 

with 




l + 2r? 
1277 

2 "V, 



12y7 / 24r/ 



.4 



1 (v-1) 2 -Z(6 V g(l;a) + l) 



12?? (2 + rj) - 2g[l; a) (rj - 1) 



i (l + 2t 7 )A 1 -|(2 + » ? ) 



2 



A 2 = g(l;ay " 1 \ (31) 

1 + 6r)g(l;a) 



and 



-i 



Z =(d-p 0P l =(1 + ^(^-1)) (32) 
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with the Carnahan-Starling expression M 

together with eqs.([l5]) and ( |l4| ) completes the model. This then reduces to the GMSA in equilibrium and can be 
seen as its natural generalization to a spatially isotropic noncquilibrium state. 

V. MOLECULAR DYNAMICS SIMULATIONS 

To determine the structure of the model granular system, I have performed molecular dynamics simulations of a 
small systems of 108 and 500 particles in 3 dimensions governed by the steady-state dynamics described by eq.([l6|) and 
subject to periodic boundary conditions. The density was taken to be n* = 0.5: high enough that finite density effects 
are important but low enough that the Enskog approximation is expected to be valid. The choice to simulate the 
steady-state dynamics, rather than to simulate the "real" dynamics of the cooling system, was made on the basis that 
the systems cool very rapidly so that the time scales involved in the simulations become very large, the velocities and 
energies very small and numerical inaccuracies due to round-off error are a significant problem. This could be dealt 
with by periodically rescaling the velocities (i.e., redefining the time unit as in 0], p^| ) but it is far more elegant 
and efficient to directly simulate the steady-state dynamics. Furthermore, the changes needed to implement this 
starting with a code for simulating equilibrium hard spheres are minimal. One point that does require attention is the 
instability with respect to the total momentum. Even if initial conditions are chosen so that the total momentum is 
zero at the start of the simulation, round-off errors lead to the spontaneous appearance of a non-zero total momentum 
which then quickly goes due to the instability. This effect is, however, benign and is easily suppressed by calculating 
the total momentum during each propagation step and subtracting (1/N) of its value from the momentum of each 
particle. 

The starting point for the simulations was an equilibrium configuration of velocities and positions. The value of 
the thermostat constant, w, was set arbitrarily. For each value of a, the simulations were "equilibrated" for a total of 
3 • f 6 collisions and then statistics gathered over a second series of 3 • 10 6 collisions . To obtain steady-state averages 
of one-body properties, quantities were time-averaged over periods of 10 4 collisions throughout the simulations; these 
samples were then treated as statistically independent estimates and their average and standard error computed [^9| . 
The errors in all quantities reported below are found to be small, less than 1%. The determination of collisional 
effects, such as in the pressure and the pdf at contact, is somewhat different. For any collisional quantity of the form 

A = £ A(si, (34) 

i<j 

the ergodic assumption gives 

r 

T 



(A) = i / dt A(t) 







1 f T 



where the first two lines integrate the total time dependence of the function A. The third line follows from a change 
of variable in the delta- function Pij(p l i j) is the relative momentum of the colliding pair immediately before (after) the 
collision and is the time at which the pair (i,j) collide (which could be imaginary or outside the range of integration 
indicating in either case that they do not collide). The last expression obviously reduces to a sum over collisions: 



= ^ V r^ + r^j (36) 




where 7 represents the colliding pair and which is the form used to evaluate the collisional part of the pdf at contact, 
jV( ^7 1) 47Tff(f) = O^Kj 5 ilij - !))• Using |q 7 • P 7 | = a |q 7 • p 7 |, this becomes 
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|q 7 -Pj\ 



1 



— 6(-q 7 -p 7 ) 



(37) 



while for the pressure one finds 



(1 + g) 
2T 



i q f ■ p 



'7 



6 (-q 7 ■ p 7 ) 



(38) 



and in both cases, the step function indicates that the expression is evaluated with the pre-collisional momenta. 
Finally, we present below determinations of the pdf for finite separations. These are determined in the obvious way 
by looping over all pairs of atoms and creating a histogram of the separations. The bin size used was 0.025 (hard 
sphere diameters) and these were compiled every 10 4 collisions and the results averaged to obtain the final histogram. 

Figure |l| shows the pdf at contact as determined from the simulations and from the collisional boundary condition. 
For a > 0.6, the agreement is seen to be good but it becomes worse values below this: furthermore, there appears to 
be a strong number dependence to the results with the larger system diverging more rapidly from the prediction. In 
both cases, simulations are only possible for a above some threshold: below this, the simulation code fails due to the 
time between collisions becoming smaller and smaller until the machine precision is reached. A detailed analysis of 
the sequence of collisions shows that this is due to a small number of particles with virtually no momentum relative 
to one another colliding over and over again - in other words, this is the phenomena of elastic collapse described by 
McNamara and Young p0]| . The threshold for the collapse is in the range 0.3 < a < 0.4 for the 500 particle system 
and 0.2 < a < 0.3 for the 108 particle system. In both cases, prior to the collapse, the value of the pdf at contact is 
several times the highest values shown in figure [|. 

Similar behavior is seen in the pair distribution function. For reference, Fig. |^ shows the pdf at equilibrium as 
determined from the 500 atom simulations and from the model: the agreement is seen to be excellent as is also the 
case for the data coming from the 108 atom system. Figures || and |J show the nonequilibrium part of the pdf (i.e., 
g(r) —g eq (r)) for the 108 atom system as determined by simulation and by from the nonequilibrium GMSA for a = 0.7 
and 0.5 respectively and in both cases, the extended MSA is seen to give a good quantitative description of all features 
of the nonequilibrium structure. (In fact, since the MD results are, by their nature, binned, the theoretical curve is 
obtained by integrating the model pdf over bins of the same size and position as used in the simulations.) Figures 
^ and ^ show that, not surprisingly, the agreement is less good for the 500 atom system, particularly at the smaller 
value of a. 

The disparity between the results for the two systems is much greater than one finds in equilibrium and suggests 
that a qualitative difference between them. One obvious possible source of such a difference is the hydrodynamic 
instability discussed above. Using the values for the transport coefficients given by |Tt|] , one finds the critical size 
curve shown in Fig. |^ which indicates that the 108 atom system is always stable but that the 500 atom system 
becomes unstable around a < 0.6 however, knowledge of the transport coefficients at small a is only approximate so 
these numbers may only be indicative of the position of the instability. Nevertheless, the importance of the instability 
in the larger system is easily confirmed and Fig. §L showing the velocities in one direction versus the positions along 
another as taken from a snapshot of the 500 atom system with a = 0.5, shows a spontaneously formed shearing profile. 
In a larger system, this would manifest itself in the form of vortices. Further evidence of the instability can be found 
in the kinetic contributions to the pressure tensor where, beginning at a = 0.7 in the 500 atom system, an oscillation 
develops whereby a large fraction, on the order of 2/3, of the kinetic energy is concentrated in first one component 
of the pressure tensor and then another indicating that macroscopic flows are forming. The obvious interpretation 
is that around a = 0.7 the shear mode is soft or unstable. There is no evidence of such an unstable mode for any 
value of at a in the 108 atom system. This picture is thus in qualitative agreement with the predictions based on the 
calculations described above and has recently been observed in other studies [Q . 

We can suppress the instability in a crude way by periodically adjusting the velocities of the atoms. In these 
"constrained" simulations we interpose correction whereby after every 100 collisions, we calculate the amplitude of 
the longest wavelength Fourier modes of the systems (i.e., A; = jj Yli=i c « cos kz • q^ for ki = (2w/L) x, etc.) and we 
then subtract the mode from each atom's velocity ( Ci — > c, — A; cosk; • q,). This is a crude procedure in that 
the amplitudes of the modes are only approximately set to zero and it also has the effect of removing kinetic energy 
from the system (which is, however, masked by the input of kinetic energy coming from the equations of motion). A 
more elegant procedure could be devised based on standard nonequilibrium molecular dynamics techniques such as 
Gauss' principle of least constraint |3l]| , but as the present purpose is only to control the unstable mode, the crude 
method was deemed sufficient. As shown in Fig. ^, the result is to give better agreement in the measured value of 
the pdf at contact between the two systems while having relatively little effect in the smaller system except at the 
highest values of alpha. Figure shows that the pdf as determined from the constrained simulation of the larger 
system is in considerably better agreement with the model. 



9 



VI. DISCUSSION 



The main purpose of this paper has been to show that the GMSA can be extended to nonequilibrium systems by 
replacing the equilibrium input required by the GMSA with accessible nonequilibrium information coming from the 
collisional boundary condition and, incidently, to elucidate the atomic-scale structure of the HCS of granular fluids. 
In order to compare the predicted values of the pdf at contact and the pressure with the results of simulations, the 
equations of motion of the dissipative hard sphere system were mapped onto those that describe a steady state thus 
allowing us to use standard methods of steady state simulation such as the replacement of ensemble averages by time 
averages. The comparisons with molecular dynamics simulations also show that the pair distribution function at 
contact can be used as a signal of the onset of elastic collapse - its value steadily diverges from the predicted value as 
the elastic collapse threshold is approached and its value in the simulations that feature the collapse is very large. 

The significant differences observed between the 108 and 500 atom systems were seen to be largely due to the soft 
hydrodynamic modes present in the larger system. Nevertheless, even when these are accounted for, there remains 
a significant deviation of the simulation results from the various predictions of the Enskog theory. It is tempting 
to conclude that this is due to a poor estimation of the quantity \o, for which we use the equilibrium value, but 
examination of the pressure shows that the answer cannot be this simple: the pdf at contact would require a larger 
value for xo that increases with a in order to be in agreement with the simulations whereas Fig. [ll] shows that pressure 
would require a smaller value that decreased with a (in both the unconstrained and constrained simulations). We 
thus conclude that the deviations are due to the Enskog approximation itself and could probably be described via 
mode coupling. Nevertheless, that xo should depend on a is intuitively clear: atoms moving slowly away from a 
collision will be more likely to be knocked, by a third atom, into a second collision with one another leading to such 
a dependence. 

As mentioned in the Introduction, one phenomcnological approach to the description of nonequilibrium structure 
is that of Hess and coworkers H , • In its simplest form, this reduces to a relaxation model for the nonequilibrium 
contribution to the pdf: 

J^(r, t) + v(r, t) ■ V<?(r, t) = t~ x ( 5 (r, t) - g (r)) (39) 

where r is a relaxation time and go (r) is taken to be the equilibrium pdf. It is clear that for a homogeneous steady state 
with no flow this gives the trivial result that g(r,t) — go(r). However, this model is intended only as a simplification 
of a more complex model given by 

J^(r, t) + v(r, t) ■ Vg(r, t) + DV ■ (g (r, t)V (ff(r, t)/g Q {v))) = (40) 

where D is a kind of diffusion constant. Again, for a homogeneous steady state with no flow, only the last term 
survives so that this reduces to 

1 d 2 d g(r) dg (r) d g(r) 

9o( r ) — ^- r -} rr + — j — rr = (41) 

r^ dr dr go{r) dr dr go{r) 



which has the solution 



f(r) = / dr—-- (42) 

where we have used the boundary condition linv^oo g{r) = 1. Although this has an interesting structure, the function 
f(r) is positive-definite so that for HCS the difference g(r)—go(r) is also which precludes a description of the oscillatory 
nature of the differences found in the simulation. 

Another approach is the theory of Eu et al [JL2J, Jl3| , also mentioned in the introduction. This theory is based on 
an ansatz for the N-body nonequilibrium distribution and takes the form of an integral equation: 

lny(qi,q 2 ;i) = n J dq 3 f NE fa, q 3 ; t) y(ct_, q 3 ; t) 

x [yfa, q 3 ; *) (i + Jne fa, qs;*)) - 1] (43) 

where 
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9 (qi , q25 t) = exp [-V NE (qi , q2! *)] V (qi, q2; 
/jvu (qi,q3;*) = exp [-V NE (qi,q2;*)] - 1 ( 44 ) 

and where the nonequilibrium potential Vn e (qi , Q2 ; t) is a sum of the equilibrium potential and terms related to the 
moments of the velocity. This formulation has the desirable property that, if linearized in the density by replacing 
lnj/(qi, q2; t) — J/(qi,q2ji) — 1) it reduces to the PY approximation in equilibrium. It is difficult, however, to see 
how to incorporate the exact requirement of the information coming from the collisional boundary condition and it is 
therefore an open question whether this can give a model comparable to the nonequilibrium GMSA described above. 
Indeed, in the derivation of this model, the authors explicitly neglect velocity correlations of the kind used here to 
control the model of the structure. Nevertheless, the solution of this model to allow for such a comparison would be 
of some interest. 

It is also appropriate to comment on the applicability of the Enskog-level description of the system. The density 
chosen for the simulations is one at which the Enskog description of elastic hard spheres is very good with most 
quantities being accurately predicted to within a few percent. Furthermore, neither the calculation of the pressure 
nor the pdf at contact requires explicit knowledge of the one-body distribution at contact which is good, because no 
exact solution to the Enskog equation for HCS exists. None the less, the results described above show that there are 
systematic deviations from the Enskog predictions at all values of the coefficient of restitution. Evidence has also 
been given that these are at least partly due to the presence of soft modes within the system. The conclusion is 
therefore mixed: while the Enskog description seems to be qualitatively accurate, all but the very smallest systems 
will contain soft or unstable modes that result in significant deviations from it. This, more than the separation of 



time scales discussed in ref. 32 would seem to be the greatest obstacle to using a Boltzmann/Enskog description 
of the one body function or of using a hydrodynamic description of the macroscopic state. Indeed, the Enskog and 
hydrodynamic descriptions are successful in predicting fairly well the location of the hydrodynamic instability in the 
vorticity. A better test of these issues would be to study a related system, such as a sheared granular fluid, which 
may be more stable. 

The nonequilibrium GMSA described here works surprisingly well. Encouraging results have also been found when 
this model was applied to simple sheared fluids Q and a systematic study of this system is in progress. It is also of 
interest to try to improve on the use of the compressibility equation in the nonequilibrium model: one substitute would 
be information coming from kinetic theory such as the asymptotic behavior of the pdf for which various approaches 
are possible. 
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APPENDIX A: ORIGIN OF THE COLLISIONAL BOUNDARY CONDITION 



1. From BBGKY hierarchy 



It is instructive to derive the collisional boundary condition from the BBGKY hierarchy, the first equation of which 
appears above as eq.(0). The second BBGKY equation is 



d d c^ — . . 

■777 +Pl ' T. hp2 ' T. T-(xx,X 2 ) 

at dqi dq 2 



f {2) {x u x 2 -t) 



where the collision operator is 



dx 3 \T_(xi,x 3 ) + T_(x 2 ,x 3 )] f ( > 3 ')(x 1 ,X2,x 3 ;t) 



T-(x 1 ,x 2 ) = cr / da 6 (pi2 • a) (pi 2 • a) 5 (q 12 - a a) a b 12 - S (q i2 + a a) 
Jq 1 

= J da S (a - a ) 6 (P12 ■ <?) (P12 • °) 5 (q i2 - a) aT 2 ^ 1 - 5 (q 12 + a) 
= S {q 12 - a ) (P12 • qi2) 6 (P12 ■ qi 2 ) a^ 2 b^ + 6 (-pia • qi 2 ) 



= S (qi2 - er ) (p i2 • qi2) 



a~ 2 K} + 1 



© (-P12 • qi2) 



(Al) 
(A2) 

(A3) 
(A4) 

(A5) 
(A6) 



We now observe that no matter what the state, the atoms cannot overlap so we must be able to write the distribution 

as 



fW( Xl ,x 2 ;t) = 9 (q 12 - a ) f^{x x ,x 2 ;t) 



12)/ 



so that 



Pi ' ^ h P2 ' 

oqi aq 2 



e(qi 2 -a )f (2) (xi,x 2 ;t) 



= P12 • qi 2 <5 (qi2 ~ o- ) r 2 '(xi, x 2 ] t) + 6 {q 12 - a ) 



Pi ' ^ 1" P2 • 

dqi dq 2 



fW(x u x 2 ;t) 



(A7) 

(A8) 
(A9) 



so that there are two sources of singularities in the second BBGKY equation: one coming from the collision operator 
and one from the streaming operator. If we integrate q\ 2 over a vanishingly small interval centered at <7o, only the 
singular terms will contribute so that we conclude they must cancel independently of the regular terms. This gives 



or 



P12 • qi2<5 (gi2 - a ) f( 2 \xi,x 2 ;t) = S {q X2 - a ) (p 12 ■ q i2 ) 



S (qu - ct )6 (pi2 • qi 2 ) f {2) (xi,x 2 ;t) = S {q 12 - a ) a 2 b 12 L e (-p X2 • q i2 ) f w {xx,x 2 \t) 



-21-1 



1 



e(- P i2-qi2)/ (2) (xi,a;2;t) (A10) 



'(2)/ 



(AH) 



which is the desired result relating the post-collisional distribution, on the left, to the pre-collisional distribution on 
the right. It is thus apparent that this identity carries part of the information of the second BBGKY hierarchy. 
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2. From Conservation of probability 

It is clear that the probability to find two atoms moving towards each other with a given relative momentum, the 
left hand side, is the same as that to find two atoms in contact moving away from one another with the corresponding 
post-collision momentum: 

<5(qi2 -<7CT)6(-qi2 • P12) h{x\, x 2 , t)d 3 x 1 d 3 x 2 = S (qi 2 - ad) 6 (qi 2 ■ p' 12 ) f 2 (x' l7 x' 2 ,t)d 3 x' 1 d :i x / 2 

which simply states that the probability to find two atoms in contact moving towards each other with a given relative 
momentum, the left hand side, is the same as that to find two atoms in contact moving away from one another with 
the corresponding pre-collision momentum (the right hand side). Using 



1 

a 2 

f 2 {xi,x 2 ,t) =biif 2 (x' 1 ,x' 2 ,t) (A13) 



d 3 x x d 3 x 2 = —d 3 x[d 3 x 2 (A12) 



this can be written as 

5 (qia - ad) 6 (q i2 • p' 12 ) \ b^ 2 f 2 (x[, x 2 , t) = 5 (q i2 - ad) 9 (q i2 • p' 12 ) f 2 (x[, x 2 , t) (A14) 

or 

. The factors of a in the Jacobian arise because of the change from pre- to post-collisional momenta (p' 12 = — api 2 
gives one factor of a) and because the relation between positions q(i) after the collision and positions before the 
collision also involves the momenta giving a second factor of a. 

3. Evaluating correlations at contact 

From the identity and the assumption of molecular chaos, we have 

<5(qi2 - aa)Vf 2 (xi,x 2 ,t) = S (q i2 - ad) /i(xi; i)/i(x 2 ; t)# (qi, q2; t) (A15) 

+5 (q 12 - ad) 6 (q i2 • p 12 ) (a^ 2 b^ - lj /i(ari; i)f 1 (x 2 ; t)g {q.i, q 2 ; t) (A16) 

so that if averages over /i(xi; i)/i(x 2 ; £)#o(qi, q2! t) are denoted by (...} , then for any two-body function A = 
Y,j<i A(x i ,x j )5 (q»j ~ da) one finds 

A(d) = (A) (A17) 
= (A) + — J dx Y dx 2 5 (q i2 - ad) A^i, 2:2)6 (q i2 • P12) {oT' 2 b^ 2 1 - l) h{x\; t)fi(x 2 ; £)#o(qi, q 2 ; t) (A18) 

= (A) + — J dx Y dx 2 S (qi2 - ad) fi(xi;t)fi(x 2 ;t)g {ca,q 2 ;t) (oT 1 ^ - lj 9 (qi 2 • P12) A(xi,x 2 ) (A19) 

= (^)o + { S ( q *J ~ ( e ' Pt ^ ~ e (%J ' P«)) A ( x h x jj) ( A2 °) 

j<i 

and, in particular, if A(xi,Xj) = then 

g(a) = g (a) + ^^9o{d) (A21) 
= (A22) 

as reported in the text. 
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FIG. 2. The pdf at equilibrium (a — 1) for n*=0.5 as determined from simulation of 500 atoms (circles) and from the(GMSA) 
model (curve). 
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Fig. 3 
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FIG. 3. The nonequilibrium part of the pdf, 5g(r) — g(r;a) — g(r; 1), for n*=0.5 atoms with a = 0.7 as determined from 
simulation of 108 atoms (circles) and from the model (curve). 
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FIG. 4. Same as fig. 3 for a = 0.5 and 108 atoms. 
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Fig. 5 
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FIG. 5. Same as fig. 3 for a = 0.7 and 500 atoms. 
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Fig. 6 
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FIG. 6. Same as fig. 3 for a = 0.5 and 500 atoms. 
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FIG. 7. The critical system sizes as a function of a for n*=0.5. Systems falling below the curve are unstable. 
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Fig. 8 
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FIG. 8. A snapshot of the 500 atom system: the horizontal axis is the position along the x-axis of the simulation, the vertial 
axis shows the momentum along the z-direction. The curve is a sin-function fitted to the data. 
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Fig. 9 
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FIG. 9. The pdf at contact for n*=0.5 from the original unconstrained simulation of 108 atoms (open circles), the constrained 
simulation of 108 atoms (circles) and the constrained simulation of 500 atoms (squares) and from eq, (|l4|) (line). The lines 
between the simulation data are only a guide to the eye. 
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Fig. 10 
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FIG. 10. Same as fig. 6 for a = 0.7 and 500 atoms and showing data from the constrained simulation. 
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